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Abstract 

This paper describes and discusses the textbook. Fundamentals of Computational 
Fluid Dynamics by Lomax, Pulliam, and Zingg, which is intended for a graduate- 
level first course in computational fluid dynamics. This textbook emphasizes funda- 
mental concepts in developing, analyzing, and understanding numerical methods for 
the partial differential equations governing the physics of fluid flow. Its underlying 
philosophy is that the theory of linear algebra and the attendant eigenanalvsis of 
linear systems provides a mathematical framework to describe and unifv most nu- 
merical methods in common use in the field of fluid dynamics. Two linear model 
equations, the linear convection and diffusion equations, are used to illustrate con- 
cepts throughout. Emphasis is on the semi-discrete approach, in which the governing 
partial differential equations (PDE’s) are reduced to systems of ordinarv differential 
equations (ODE’s) through a discretization of the spatial derivatives. The ordinary 
differential equations are then reduced to ordinary difference equations (OAE’s) us- 
ing a time-marching method. This methodology, using the progression from PDE 
through ODE s to OAE’s, together with the use of the eigensystems of tridiagonal 
matrices and the theory of OAE’s, gives the book its distinctiveness and provides a 

sound basis for a deep understanding of fundamental concepts in computational fluid 
dynamics. 
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enceoterato^ P T S T^rl'i *" ^ Paper ’ include the model ec I u ations, matrix differ- 
perators Tavlor Tables, the semi-discrete approach, converting time-marching 

methods to ordinary difference equations (OAE’s), solution of linear OAE’s the A -a 

of ^he°m r d " UmeriC !i 1 Stab * h ^ V conce P ts - Although this omits a significant portion 

Most of the" 13 f C0Ve !; ed m the book ’ 11 13 sufficient to illustrate the basic framework. 
•Most of the material presented is excerpted directly from the book. 


2 The Model Equations 

uleli C ° nCeptS are P resentecl in ^e context of two one-dimensional model equations 
mear convection equation and the diffusion equation. These two equations are 

scalar °and thl m - Under ® ta " din « CFD algorithms, primarily because they are linear, 
scalar and they represent phenomena of importance to the analysis of certain aspects 

of fluid dynamics. The linear convection equation can be written as- 


du du 

aF + “S =0 


(i) 

of involvmg convection " ithout considera,ion 

du d 2 u 

~dt ~ U dx^ (2) 

u here i/ i S a positive real constant. Consideration is given to both Dirichlet (specified 
“) and Neumann (specified du/dx) boundary conditions. leaned 

3 Finite-Difference Operators 

Fimte-difference operators are presented initially through linear combinations of Tav- 

erltor fnrT PanS, T!i’ P Ci " g ’ for exam P ie . ‘he three-point centered difference op- 
erator tor a second-derivative: ^ 


(S Ix u)j - ——( Uj+l _ 2 uj + 


( 3 ) 

The idea of a local order of accuracy is also defined based on the leading term in the 
Taylor series expansion* of the error. 

The discrete Taylor series expansion is defined as = Uj + (kAx)(^) j + L(kAx )- + 
••• + j(l:Ax) ,l (|J). + • •• for a general increment k ‘ 
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notation is used: 


B P (M : a, b , c) = 


b c 
a b 


a b c 

a b 


( 8 ) 


M 


With this notation, the matrix difference operator for a first derivative with periodic 
boundary conditions on an i\/-point grid can be written as 


1 




(9) 


Notice that there is no boundary condition vector (be) since this information is inte- 
rior to the matrix itself. 


3.2 Taylor Tables 

A Ta\ lor table is a simple and convenient way of forming linear combinations of 
Tavlor series on a term by term basis. This enables one to derive finite-difference 
expressions for specified derivatives with a specified stencil. An example is given 
below. The table is constructed so that some of the algebra is simplified. At the 
top of the table we see an expression with a question mark. This represents one of 
the questions that a study of this table can answer; namely, what is the local error 
caused by the use of this approximation 7 Notice that all of the terms in the equation 
appear in a column at the left of the table (although, in this case, Az 2 has been 
multiplied into each term which simplifies the terms to be put into the table). Then 
notice that at the head of each column there appears the common factor that occurs 
in the expansion of each term about the point j, that is, 


Az* 



k - 0, 1,2, 


The columns to the right of the leftmost one, under the headings, make up the Taylor 
table. Each entry is the coefficient of the term at the top of the corresponding column 
in the Taylor series expansion of the term to the left of the corresponding row. For 
example, the last row in the table corresponds to the Taylor series expansion of cv,j+\: 

' du\ 


1 


cu j+l = CUj -t-c- (1) • jAx ■ +C- (l) 2 - ^Az 2 

c-( l) 4 - — Az 4 
v ' 24 


+c- (l) 3 • jr Az 3 




( 10 ) 
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4 The Semi-Discrete Approach 


equation (PDE) is m ? J | e 'i llfference approximations to a partial differential 
v r r ™ r by dlfferenClng the s P ace derivatives only, therebv con- 

tin 0 the PDE to a set of coupled ordinary differential equations (ODE ! s) We refer 

to th,s as the semi-discrete approach. For some methods it is not possible to sepa- 
rate the spatial and temporal discretization, and two of these, the Lax-Wendroff and 
. lacCormack methods, are discussed later in the book. However, we primarilv con- 
centrate on the semi-discrete approach, reducing the PDE’s to ODE’s bv discretizing 

he develo 7 7 ^7 ^ deVel ° ped theor >’ of 0DE a^tioni to aid us in 
the development of an analysis of accuracy and stability. 

„ n( W' he m . odel PDE ’ S ’ " e can approximate the space derivatives with difference 
operators and express the resulting ODE’s with a matrix formulation. This is a 

central diffe™" at '° n u th<! 0DE ' S are linear ' For example, using the 3-point 


da v _ 

(be) ( 13 ) 


with Dirichlet boundary conditions folded into the (be) vector. 

The generic matrix form of a semi-discrete approximation is 
equation 


expressed by the 


du 

di = Au ~ nt) (14) 

HifiWp hat the elements in the matrix depend upon both the PDE and the tvpe of 
t iZ n boLt eme hT th , 6 SPaCe t6rmS - The VeCt ° r is USuall - v determined 

aL Wr Sto^ C ° nd, ,0nS and r ib,y S ° UrCe termS ‘ In general > — the E ^r 

a ler-Stokes equations can be expressed in the form of Eq. 14 In such cases 

** usuatv d S ™ ^7 7 ^ ° f A depend ° n the elution ff and 

allv derived by finding the Jacobian of a flux vector. Although the equations 

are non mear, linear analysis leads to diagnostics that are surprisingly accurate when 

° f nUmerical methods as they apply to the Euler and Navier- 

nf 1,177 th 7 U *a M matrLX ' 4 iS inde P endent of « ^ d t and has a complete set 
• I [{ mde P endent eigenvectors, the generic system can be written as a set of 
independent first-order equations in the form: 


w\ = XiWi-g^t) 
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The numerical solution to a set of linear ODE’s (in which A is not a function 
of t) is entirely equivalent to the solution obtained if the equations are decoupled, 
solved in uncoupled form, and recoupled. The importance of this concept resides in 
its message that we can analyze time-marching methods by applying them to a single, 
uncoupled equation and our conclusions will apply in general. This is helpful both 
in analyzing the accuracy of numerical time-march methods and in studying their 
numerical stability. 

As a single representative ODE we choose: 


— = Xu + ae» ( 21 ) 

where A, /.i, and a are complex constants, which can be used to evaluate all manner 
of time-marching methods. In such evaluations the parameters A and f.i must be 
allowed to take the worst possible combination of values that might occur in the 
ODE eigensystem. The exact solution of the representative ODE is (for // ^ A): 

n(t) = ce M + (99) 

f.i - A v ' 

Although it is quite simple, the numerical analysis of Eq. 21 displays many of the 
fundamental properties and issues involved in the construction and studv of most 
popular time-marching methods. 


5 Time-Marching Methods for ODE’s 

5.1 Converting Time-Marching Methods to OAE’s 

As we have seen, application of a spatial discretization to a PDE produces a coupled 
system of ODE’s, which can be solved numerically using a time-marching method. 
Application of a time-marching method to an ODE produces an ordinary difference 
equation (OAE). In this section, we shall analyze and contrast explicit and implicit 
schemes to illustrate these methods. For example, consider the explicit Euler time- 
marching method, which can be written as: 


where 


^n + 1 — -f- hu n 



F(u,t ) 


(23) 


(24) 
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5.3 The X-cF Relation 

io ^ft f r of 

time. The exact solution is PProach, Eq. 14, with / independent of 

^ = c,R)"*- 1 + ... + Cm ( e ^ r S m + ... +eu ^y iu + A -, f {33) 

Where t = nh, while the explicit Euler method produces the solution 

u„-c,(,7,) x, + " • + c m ((7„)" X m H +Ci,(ca)“i„ + ES. (34) 

between and e- Since the va Jof^t, 


h 


_ 1 + Ah + ^X 2 h 2 + -U'V + • • • + ~\ n h n + ■■■ 

D n: 


(35) 

the truncated expansion <7 = I + AA is an approximation to e« valid for small AA 

m a TT “ a "d° f “ he 

for every A eigenvalue that satisfies 1 re lat on ” ^ 


o - 1 + AA + i A V + • - . + I A*A* + 0(A‘7') (36) 

where k is the order of the time-march method. We refer to the root that ha* ,h 

deta'ils^r the') aS ““ CT - r00t The property can be stated regardless of the 

he time-march method, knowing only that its leading error is 

hus the principal root is an approximation to e xi up to 0(h k ). When the time 
march, ng method incorporates the solution at time level n - 1 !r earlier (for exam!" 
he leapfrog method), add.tional o-roots arise, known as spurious roots These nlav 

U Le oTpTh ana ' ySiS ' , bU ‘ m “ St bC C °" Sidered in stabili ‘y analysis. 

matrix A and E \ ” ““n" T"' 60 " the eigerlvalues of the spatial operator 

The Ai^z't' 

solutton can also be determ.ned from the 04E solution. Exampfes of t"s 


& — 1 + A h 


(37) 
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2. \\ ill the numerical solution converge to the exact solution of the PDE as the 
grid is refined? 

The second question leads to the notion of Lax stability and Lax’s theorem, which can 
quickly lead to analysis that is quite complex for a first course. Hence we concentrate 

on the first question, that is, we approach stability primarily from the ODE point of 
view. 

This approach leads to the criterion 1 


km I < 1 for all m 


(44) 


for both principal and spurious roots, i.e., all roots must lie on or within the unit 
circle in the complex plane. Using the A -a relation for any time-marching method, 
one can plot its stability contour, which is the locus of points for which the largest |cr|> 
considering all a- roots, both principal and spurious, is equal to unity. This contour 
is plotted in the complex A h plane and divides this plane into stable and unstable 
regions. For a given ODE in the form of Eq. 14, the spectrum of eigenvalues of A 
must lie within the stability contour of the time-marching method. This requires the 
choice of an appropriate method and time step. 

For example, consider the family of one-step time-marching methods given by the 
following expression 


u n + 1 — U n + 


~ e ) U 'n + 0U n + l] 


(45) 


This family includes the explicit Euler (9 = 0), the trapezoidal (, 9 = I), and the 
implicit Euler methods (9 = 1). Its A — cr relation is 


a = 


1 + (1 — 6)Xh 


(46) 


i - exh 

Fig. 1 displays the corresponding stability contours. The explicit Euler method is 
stable for a finite portion of the left half-plane. This leads to conditional stability. 
For a given set of A-eigenvalues, there is always an upper bound on the time step 
abo\e which the method is unstable. Note that for imaginary eigenvalues the explicit 
Euler method is unconditionally unstable. Hence it is inappropriate for application 
to the linear convection equation with periodic boundary conditions and centered 
space differencing, since this leads to pure imaginary A-eigenvalues, as shown in Eq. 
20. Both the trapezoidal and implicit Euler methods are stable for the entire left 
half-plane and are thus unconditionally stable. These examples show typical stability 
characteristics for explicit and implicit methods, which must be weighed against the 
increased cost per time step associated with implicit methods. The deciding factor is 
often the stiffness of the ODE system being solved. 


* Actually the criterion is \a r 
book. 


< 1 if one includes defective systems, which are covered in the 
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